BEST CONSTANTS IN LIEB-THIRRING INEQUALITIES: A NUMERICAL 

INVESTIGATION 



ANTOINE LEVITT 



Abstract. We investigate numerically the optimal constants in Lieb-Thirring inequalities by 
studying the associated maximization problem. We use a monotonic fixed-point algorithm and 
a finite element discretization to obtain trial potentials which provide lower bounds on the 
optimal constants. We examine the one-dimensional and radial cases in detail. Our numerical 
results provide new lower bounds, insight into the behavior of the maximizers and confirm some 
existing conjectures. Based on our numerical results, we formulate a complete conjecture about 
the best constants for all possible values of the parameters. 



1. Introduction 



In this paper we study the family of Lieb-Thirring inequalities [20,21 , which state that for 
any potential V G L 7+ 2(M d ,IR) and non-negative exponent 7, 

tr((-A + V)I)<L 7jd J Vl + -\ (1) 

where V- = max(— V, 0) is the negative part of V, and (—A + V)Z. is the 7-th power of the 
negative part of the operator (—A + V), as defined by the functional calculus. In other words, 
tr((— A + V)Z.) = J2i l^il 7 i s the 7-th moment of the negative eigenvalues of —A + V. 

This inequality was originally used by Lieb and Thirring in the case 7 = 1 , d = 3 as an 
important tool to prove the stability of fermionic matter in |20] (see [17 18 for an overview). 



The generalization ([!]) to general 7 and d has since then attracted a great deal of attention (see 
for instance [IT] for a review). Of particular interest are the values of 7 and d for which this 
inequality holds, and the value of the optimal constants L 7 d- 

Despite the physical significance of the Lieb-Thirring inequality and the amount of mathe- 
matical research on the subject, many questions remain. This paper aims to investigate some 
of these numerically, something that, to our knowledge, has not been done since the work of 
Barnes in the appendix of the original paper by Lieb and Thirring [2l], in 1976. 

We now briefly review what is known about L 7i d- It is easy to see that the bound ([!]) cannot 
hold for 7 < 1/2 when d = 1 by scaling, and 7 = when d = 2 (because any arbitrarily 
small potential in two dimension creates at least one bound state). The inequality was proved 
for 7 > 1/2 in one dimension and 7 > in other dimensions in the original paper by Lieb 
and Thirring |21 . The proof in the case 7 = 1 was recently greatly simplified by Rumin |25| . 
The case 7 = 0, d > 3 requires completely different methods and was proven independently by 
Cwikel, Lieb and Rosenblum [3j[l5j[24]. The limit case 7 = 1/2, d = 1 remained unsolved for 
twenty years until it was finally settled by Weidl [27], and refined with a sharp constant by 
Hundertmark, Lieb and Thomas [9]. 

The inequality ([T| can be interpreted as a comparison of the energy of the quantum mechan- 
ical system given by its Hamiltonian — A + V, with its semiclassical approximation that only 
involves integrals of V. The semiclassical regime is obtained by letting the Planck constant h 
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tend to zero in the Hamiltonian —hA + V. In this paper, we have scaled h away from inequality 
for convenience, and the corresponding limit is a large (or extended) V. The eigenfunc- 
tions become localized in the region of the phase space defined by p 2 + V(x) < and explicit 
computations are possible in the limit. More precisely, using the Weyl asymptotics, one can 
prove 

r tr((-A + jiV)!) _ r 

11111 ~d ~~ ^y,d,sc, 

where 

T - o-d„-d/2 r (7 + 1) 

7 ' d ' sc " r( 7 + f + i)- 

Therefore, denoting by L 7i <j the optimal constant in ([!]), we have L 7i ^ > L 7i d iS c, and we set 

R l4 = > I- 

The value of Rjd describes by how much the quantum mechanical energy can exceed its 
semiclassical counterpart. This paper is an investigation of R^ }d . It aims at providing insight 
as to its behavior by trying to solve numerically the variational problem 

R l4 = j^- sup tr((-A + V)Z), (P J>d ) 

where we impose the condition J V y+ 2 = 1 to eliminate the scaling invariance of the Lieb- 
Thirring inequality. 

Several important facts about Ry,d are known. A scaling argument by Aizenman and Lieb [TJ 
shows that R l4 is decreasing with respect to 7. Laptev and Weidl [12] showed that R Jid = 1 
for 7 > 3/2. Helffer and Robert |8] proved that R l4 > 1 for 7 < 1 by expanding tr(— A + V)Z 
for the harmonic potential V{x) = 1 — |x| 2 in the semiclassical limit. From these results, we can 
deduce that for each dimension d there exists a critical 7 Cj d such that 

<c,d, 



I R ltd > 1 when 7 < 
\Ry,d = 1 when 7 > 7 Cjd , 



with 

3 

1 < lc,d < 2- 



A trial potential of 21 provides a lower bound on Rjd an d therefore 7 C . In the appendix of 



this paper, Barnes solved numerically the restricted problem 

L 



Ry,d,l = Y~ SU P ( 2 ) 

' 7 ,d,sc Jy7+| =1 



where Ai is the lowest eigenvalue of —A + V. This is equivalent to restricting (Pj t d) to the 
potentials V such that —A + V only has one negative eigenvalue. The solution V^,d y i of this 
problem is negative, radial and only has bound state, i.e.— A + Vy idj i only has one negative 
eigenvalue. The corresponding Rj,d,i is decreasing, and intersects 1 at a critical 7^7. In low 
dimensions, 7 c ,i,i = |, 7c,2,l ~ 1-165, J c ,3,l ~ 0.863 (our numerical results agree with these 
values). Therefore, 7 Cj i = | and 7 Cj 2 > 1.165, but nothing can be said about 7 c .d for d > 3. A 
famous conjecture of Lieb and Thirring states that R~ d, = 1, ie.7 C) 3 = 1. This would imply 
improved bounds on the energy of a system of fermions, and is of great importance to the 
Thomas-Fermi theory |16| . 

Better upper bounds on R™ d have also been derived recently. For instance, it is proven that 
R-y,d < ~ 3.63 for 7 > 1, R*y t d < ^5 ~ 1-82 for 7 > 1 (see |6j), by generalizing the inequality 
to matrix- valued potentials. 



Despite these advances on upper bounds, not much is known about lower bounds for FL, ^, 
except for the one-bound state potential Vjdi of [2l] and the asymptotic result in the semi- 
classical limit of pi. In this paper, we attempt to bridge the gap between these two results by 
looking numerically for maximizers of the variational problem (Py^). To our knowledge, this 
is the first numerical study of the Lieb-Thirring inequalities since the work of Barnes |21j (and 
unpublished recent work by Arnold and Dolbeault W in ID). 

The method we use is a finite element discretization of a fixed point algorithm. We describe 
this algorithm, and specialize it to the case of radial potentials. Then we discretize it, and use 
it to obtain numerical results. The critical points we obtain are the natural generalization of 
the potential with one bound state obtained by Barnes. They serve as a partial bridge between 
this potential and asymptotic results, and yield new lower bounds on P 7) d. 

2. The optimization algorithm 
2.1. Optimization scheme. Let us denote E(V) = tr((— A + V)Z), so that 

L-y,d= sup E(V). 

The maximization algorithm is based on the following well-known property: 
Proposition 1. For any V G L 7+ 2(M d ,R) and any j > 1, 

EiV) 1 ^ = max -tr((-A + V», 

T>0 

IM|y=l 

/ \ 1/7' 

where ||r||y = ( tr (j T l 7 ') J ?s the Schatten norm of the operator r with exponent 7' = 7 zy, 
the Holder conjugate ofj. 
Proof. We have 

- tr((-A + V) t) < tr((-A + V)_ r) 

<||(-A + V0-|| 7 ||r|| y 

< trd-A + V)!) 1 ^ 

< P(y) 1/7 , 

and the equality is achieved when 

r = K T (-A + V)T\ (3) 



where K T is a normalization constant, chosen to ensure that ||r|| / = 1. 



□ 



From the proposition, we see that when 7 > 1, 

J i4 



(L 7|d ) 1/7 = sup -tr((-A + y)r) 



r>0,||r|[y=l 

What we gain from this formulation is that we can explicitely maximize — tr((— A + V)r) 
with respect to V and r separately. Indeed, for a fixed V, the maximum with respect to r is 
given by (|3]), and for a fixed r, the maximum with respect to V is given by 



V(x) = -KvTix.x)^^ 



where Ky is a normalization parameter chosen to ensure that fV~ 1+ 2 = 1, and r(x,y) is the 
integral kernel of r. 

This suggests the following maximization scheme: Given an approximate maximum V n , we 
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set T n = K T (— A + V n )_ and V n +i(x) = —KyT n (x,x)~' + 2- 1 l and iterate. Explicitely: 
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Algorithm 1. Maximization algorithm 

(1) Compute the negative eigenvalues Aj and corresponding eigenvectors ipi of —A + V n 

(2) Set Pn = EiC-Ai) 7-1 ^ 2 





i 


-1 




(3) Compute K n = 


7+1-1 
Pn 










7+1 





7+4-1 



(4) Set V n +i = -K n p n 
By construction, this algorithm ensures that 

£(K+i) 1/7 = -tr((-A + V n+ l)T n+1 ) 

> -tr((-A + y n )r n+ i) 

> -tr((-A + K)r n ) 
= E{V n ) 1 ^ 

Therefore, the sequence E(V n ) is non-decreasing, and since it is bounded from above, it 
converges. In general though, V n may not converge. We give examples in Section [4] where, 
because of the translation invariance of the functional, V n splits into two bumps that separate 
from each each other. Even when this behavior is forbidden, for instance by imposing a finite 
domain as we do in numerical computations, a rigorous convergence analysis of the algorithm 
is still an open problem. 

Note that, even if Algorithm 1 was derived from Proposition [T] for 7 > 1, it remains a 
reasonable algorithm when 7 < 1. In this case, though, the monotonicity of E(V n ) is not 
guaranteed, and indeed was found to be false in numerical computations. Nevertheless, we did 
not notice any specific convergence issue for 7 < 1. 

Algorithm 1 can also be seen as a fixed-point scheme for the critical points of the maximization 



problem. Indeed, the Euler-Lagrange equations for (P~ /j( i) are 

-Ky [(-A + V)l 



V(x) 



~ l (x,x) 



7+f 



(4) 



This is a self-consistent set of equations similar to systems such as the Hartree-Fock equations of 
quantum chemistry [19] . Our algorithm is similar in spirit to the Roothaan method [23] . In our 
case, at least for 7 > 1, the scheme monotonically increases the objective function. Therefore, 
the oscillatory behavior often seen in the Hartree-Fock model [2 14 cannot occur here. Even for 
7 < 1, we did not see any such oscillations numerically, and always observed linear convergence 
(i.e-llV^ — Voo\\ < v n for some v, < v < 1), or slow separation of bumps, as will be discussed 
in Section 01 

This algorithm is also used in the context of the Lieb-Thirring inequalities by Arnold and 
Dolbeault in ID [3]. 

2.2. Radial algorithm. Most of our multidimensional computations were done in a radial 
setting. To see why this is appropriate, consider the above iteration for d > 2, when V n is 
radial. Then, the Laplacian splits into A r + ^LtA#, where A r is the radial Laplacian, and 
Ag the Laplace-Beltrami operator on S d ~ l . The eigenvectors of Ag are explicitly known to 
be the functions whose orthoradial part are the spherical harmonics, labeled by the integers £ 
and m. Since V n is radial, —A r + V n commutes with Ag and can therefore be diagonalized in 
the same basis (separation of variables). We write these eigenvectors of the form ipi/ tT n( x ) = 
4>i,i(\z\)Je,m( x /\ x \)i where J^ m is the spherical harmonics of degree I and order m. The radial 
parts (pi £ satisfy the equation 



1 



„d-l. 



+ 



£(£ + d 



2) 



4>i,£ + V4>i,i = \,e<j>i,e, 



(5) 



and each Aj a has multiplicity 



h{d,£) 



d + i-l 




(see (26].) 

Therefore, we can obtain all the negative eigenvectors and eigenvalues by solving only ([5]). 
Furthermore, since the lowest eigenvalue decreases as i increases, we can iterate over t and stop 
whenever the lowest eigenvalue becomes positive. Next, we compute 

Pn(x) =^2^2^2(-k,e)' r '~ 1 AAm(x) 2 

l i m 

= 5>(^)E(-W~ 1 <M0 2 



where the sum only ranges over all negative eigenvalues. This is again radial, and therefore so 
is V n +i- To summarize: 

Algorithm 2. Radial maximization algorithm 

(1) Compute the negative eigenvectors <fii t £ and eigenvalues A^ of ([5J by increasing £ until 
the lowest eigenvalue Ao i becomes positive 

(2) Set p n = ZtEiHdJX-X^r- 1 ^ 



(3) Compute K n 



1 


-1 




Pn 




J P n 




T+l 


\ J 



7+f-l 



(4) Set K+i = ~K N p n 

Note that this algorithm is not an approximation: the iterates generated by this algorithm 
coincide with the ones from Algorithm 1 when the initial guess Vo is radial. 

3. Discretization 

3.1. Galerkin basis and weak formulation. To discretize this algorithm, we use a Galerkin 
basis on which we expand V and the eigenvectors. This discretization is variational at two 
levels. First, the choice of eigenvectors is restricted to a finite-dimensional subspace. By the 
min-max principle, this leads to overestimation of the eigenvalues and therefore underestimation 
of E(V) for a given V. Second, the choice of potentials V is also restricted, decreasing the value 
of max + d E(V). 



Our choice of discretization, rather than for instance finite differences, has the advantage 
that numerical examples can provide lower bounds on the best constants with machine 

precision as the only source of errors. A disadvantage is that the algorithm we use involves 
taking powers of functions. There is no exact way to do that in a Galerkin basis, and we must 
use an approximation, which causes a loss of accuracy in the fixed-point algorithm. 

For the non-radial ID and 2D cases, we simply use standard finite elements. The radial case 
is less standard, and we must derive a weak formulation of JHJ. There are several possibilities 

here. The simplest is to use a change of variable <p(r) 
back to the more standard form 



r 2 



to transform the equation 



d-l 
2 



)(* + ¥) ^ n \ 



X(p. 



We obtain a weak formulation by multiplying by a test function u and integrating: 



6V + 



+ v\4)u 



A / (f>u, 



(6) 
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Expanding the function <j> on a Galerkin basis <f>(r) = J2iXiXi( r )i this problem transforms 
into the generalized matrix eigenvalue problem 



Ax = XMx, 



(7) 



where 



Aij = I XiX'j + 



Mij = / XiXj, 



-5 + ^ XiXj, 



(8) 
(9) 



are the stiffness and mass matrices. 

When V is expanded on the same basis, we can compute the matrix elements, solve the 
eigenvalue problem ([7]), and transform back to <fi. However, for d = 2,1 = 0, ip behaves like y/r 
at 0, and the singularity of the derivative prevents an accurate discretization. 

Another possibility is to obtain a weak formulation of ^ directly. We have to remember that 
since the <fi functions are only radial parts of a d-dimensional function, the L 2 inner product 
between two functions <j>\ and cj>2 is / <i>i^>2f d ~ l . This has to be taken into account in the weak 
formulation in order to obtain a self-adjoint equation. Multiplying ^ by r d ~ 1 u, where u is a 
test function, and integrating by parts, we obtain: 



„d-i 



»'«' + 



1(1 + d- 2) 



+ V\r 



d-l 



„d-l 



m. 



(10) 



Then, as with ([6]), we can transform this into a matrix equation and solve it. The disadvantage 
is that the integrations required for the assembly step are more involved, and therefore we only 
use it for the case d = 2, / = where the other method fails. 



3.2. Finite elements. We use a Finite Element basis of piecewise linear functions on a grid 
of [0,L]. The grid we use is a nonuniform grid of N points, with more points around 0, to 
accommodate for the singularity of pi). L is chosen large enough so that all the eigenvectors 
associated to negative eigenvalues of —A + V can be represented accurately in the basis. But 
if ip is an eigenvector with negative eigenvalue A, then, ip decays as exp(— \J — Xr). This shows 
that the discretization is problematic for eigenvalues close to zero, as we shall see in Section [4] 
For our piecewise linear basis functions, it is easy to compute the matrix elements ( 8|9 ) when 
V n is expanded on this same basis Xi- We then solve the eigenproblem Q, and obtain the 

eigenvectors. To expand p 1+ i~ 1 on the basis, we simply chose the expansion that is exact on 
the grid points. Then the normalization can be performed exactly, and the iteration is carried 
out. 

We present convergence results in Figure [TJ As expected from piecewise linear basis functions, 
we obtain 0(1/ N) convergence of the eigenfunctions 4> in H 1 norm, and 0(1/ N 2 ) convergence 
of the eigenvalues A (and therefore of E(V)). For a given stepsize, the discretization error is 
exponential with respect to L, with a decay rate equal to the decay rate of the eigenfunctions, 
i.e.\/—\. Although we used a uniform grid in one dimension in Figure [IJ similar results were 
checked to hold for our non- uniform grid in d dimensions, using the weak formulation ([6]), except 
for d = 2, 1 = 0, where the convergence was slower and the weak formulation (10) had to be 
used to get the same convergence rates. 
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Figure 1. Convergence analysis with respect to N and L. The figure on the 
left illustrates the convergence of the eigenfunctions and eigenvalues with respect 
to the number of grid points N, with a fixed L = 40, large enough to cause 
a negligible error. The slopes, found by linear regression, are —1.043 for the 
eigenvectors and —1.999 for the eigenvalues, close to their theoretical values of — 1 
and —2. The figure on the right illustrates the convergence of the eigenfunctions 
and eigenvalues with respect to the domain size, with constant stepsize h = 
5 x 10~ 4 . The slope for the eigenvectors is —0.5285 « \/— A = —0.5283. When 
the domain size gets large enough, the error due to h dominates, and causes a 
plateau. 



3.3. Diagonalization. The most computationally challenging step in our algorithm is the prob- 
lem of computing all the negative eigenvalues and associated eigenvectors of a generalized eigen- 
value problem Ax = Mx, where A and M are large sparse symmetric matrices. The spectrum 
of A consists of n negative eigenvalues, and a large number of positive eigenvalues, which can 
be seen as perturbations of the spectrum of the discrete Laplacian. 

In order to compute the n negative eigenvalues, where n is unknown, we compute the k lowest 
eigenvalues, check if the k-th. one is positive, and repeat the process with a larger k if not. The 
computation of the k first eigenvalues can be done by standard packages such as ARPACK |13| . 
However, standard algorithms such as Arnoldi iteration are not suited for the computation of 
inner eigenvalues, and one has to solve for the largest eigenvalues of the shifted and inverted 
problem (^4 — a)~ 1 x = (A — a)~ l x instead to ensure reasonably fast convergence. This requires 
an adequate shift (one that is close to the bottom of the spectrum of ^4) . Even with this shift- 
and-invert strategy, the group of eigenvalues one needs to locate is not well-separated from the 
rest of the spectrum, and becomes less and less so as L increases. This leads to slow convergence 
for large L. 

3.4. Implementation. We implemented the algorithm in Python, using the Numpy/Scipy 



libraries [10], with the ARPACK eigenvalue solver 13 . 



4. Numerical results 



We used the algorithm described above to compute critical points of the variational problem 
(Py,d), that is, a solution to the self-consistent equation Q. Our strategy was to find a critical 



point by running the algorithm on a suitable initial potential Vq (for instance, a Gaussian of 
specified width). Once convergence is achieved for a specific 7, we can run the algorithm for 
7 + A7, using as initial potential the one found for 7. We have been unable to analytically prove 
the correctness of this method, for instance by checking the conditions of the implicit function 
theorem. However, we found it reliable enough for our purpose, as long as A7 was chosen small 
enough. 

7 



4.1. The ID case. In one dimension we simply used a grid of size [— L, L] with Dirichlet 
boundary conditions. We reproduced the potential with one bound state Vy t x t i of [21 by using 
for Vo a Gaussian of relatively small variance (see Figure|2]). In this case, the algorithm converges 
linearly (£.e.|| V^+i — V n \\ ~ u n , for some convergence rate v < 1), and very quickly (about twenty 
iterations to achieve machine precision). 




Figure 2. Linear convergence from a Gaussian initial data towards V~ d,\- This 
plot is for 7 = 1.2. 



Contrary to what we observe in higher dimensions, initializing the algorithm with a Gaussian 
of larger variance did not make the algorithm converge to other critical points. Instead, it leads 
to a slow divergence where "bumps" separate, each bump corresponding to the potential of 
Vy t i t i (see Figure [3]). 




n 



Figure 3. Divergence from a sum of Gaussian bumps at 7 = 1.2. Initializing 
to a single Gaussian of large width yields similar results. The two bumps sepa- 
rate slowly, until the finiteness of L forces convergence. Although not displayed 
here, the asymptotic logarithmic divergence (11) can be checked graphically, e.g. 
the spacing between Vxoo and ^1000 is the same than between V1000 an d ^10000 • 
The asymptotic slope of the log- log convergence plot is —1, which fits with the 
heuristic arguments given. 



This effect occured regardless of the value of 7 in the range 7 £ (1/2,3/2). The divergence 
appears to be logarithmic. More precisely, a numerical fit showed the asymptotic relationship 
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for the distance L n between the two bumps 

Ln ~ 



1 



= logn, (11) 

where A is the unique negative eigenvalue of —A + ^7,1,1- 

This strikingly simple relationship can heuristically be understood by the fact that, because 
the eigenfunctions <j>\ and 4>2 corresponding to the two bumps have exponential decay with 
decay constant y/—X, their interaction is of order exp(— \f— XL n ). Due to cancellations, this 
interaction leads to a correction of order exp(— 2yf— XL n ) in V n +i, and we have the approximate 
relationship L n+ \ ~ L n + K exp(— 2i/ — XL n ), which yields (11). A rigorous explanation of this 
is an interesting question. 

Based on these results, we conjecture that V^^i is, up to translation, the only maximizer of 
the functional. This would imply that 

7-1/2 



R 



7,1 



7-1/2 
7 + 1/2. 



for 7 < |, in agreement with the original conjecture of Lieb and Thirring |2lj. 

4.2. The 2D case. Due to the high cost of accurately solving the maximization problem in 
more than one dimension, we only obtained partial results in the non-radial 2D case. The results 
indicate that the algorithm either converges to radial critical points, or form slowly separating 
bumps, as in one dimension. We have been unable to find any example of non-radial critical 
points, although this does not mean that they do not exist. An example of separation in bumps 
is provided in Figure [4j 
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Figure 4. Separation in bumps of 2D non-radial initial data at 7 = 1.2, com- 
puted from Algorithm 1 with standard FEM. The initial data is taken to be a 
Gaussian in r multiplied by an angular factor (1 + cos(40)). 



In the radial case, we followed branches of critical points of (Pw) using the continuation 
method on 7 we described in Section [4} We found this branch continuation procedure robust. 
By varying the shape of the initial data, and in particular its width, we were able to obtain 
different branches of critical points. Generally speaking, as the width of the initial potential 
increases, the number of negative eigenvalues of —A + V increases. 

We display the energy of these branches as a function of 7 in Figure [5j 

9 
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7 

Figure 5. R(V) as a function of 7 for d = 2 using the branch continuation 
procedure in the radial setting (Algorithm 2). The branches are labeled by the 
number of bound states they have. Based on these results, the maximizer seems 
to be Vry t d t i, up to 7 « 1.16. Above this, no branches were above the semiclassical 
regime R = 1. 



k 


1 


4 


6 


8 


11 


13 


17 


28 


54 


81 


139 




1.165 


1.150 


1.141 


1.135 


1.126 


1.124 


1.124 


1.119 


1.110 


1.105 


1.100 



Table 1. Values ^ c ,2,k of 7 at which some branches with k bound states cross 
the threshold R = 1. 



First, we see that as the number of eigenvalues increases, R(V-y) tends to 1, the semiclassical 
limit. For a given 7, the Lieb-Thirring constant will be given by the supremum of R(Vy) over 
such curves. From the branches depicted here, we see that the supremum is always given by 
either V~ d 1 for 7 < 7* ~ 1.16, and by the semiclassical limit for 7 > 7*. All the other curves are 
below these two (although not depicted in this zoomed plot, this holds true for 0.5 < 7 < 1.5). 

It is important to keep in mind that these branches only represent critical points of the 
functional. They are generically local maxima with respect to radial perturbations, but might 
not be stable with respect to non-radial ones. For instance, Figure [6] depicts what happens 
when the radial potential with eight bound states is perturbed non-radially by a Gaussian 
multiplicative noise, with 7 = 1.1. Because the potential with one bound state has a higher 
energy, the potential splits into eight bumps. Performing the same experiment 7 larger than 
about 1.17, where the potential with one bound state has a lower energy, when the potential 
with eight bound states has a larger energy, we found that no splitting occurs, which suggests 
that the potentials are stable beyond their crossing with the one with one bound state. 

10 




Figure 6. Separation in bumps of the randomly perturbed radial potential with 
eight bound states, 7 = 1.1. 



Based on our numerical results, we conjecture that 7 C = 7* ~ 1.16. The only way this 
conjecture could be false is if some other curve is above the one with one bound state. We have 
not been able to find such a curve. 



4.3. The 3D case. Due to the high cost of computation, we were unable to get meaningful 
results in the non-radial setting, and only present our findings in the radial case. Some of the 
radial potentials we found are presented in Figure [7j with numerical data about the crossings 
of the curves in Table [21 
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Figure 7. R(V) as a function of 7 for d = 3, with the same methodology as in 
Figure [5j As 7 increases, different branches become maximizers, until 7 = 1. The 
fact that only some branches are above the threshold R = 1 when the number of 
bound states increases results from the highly oscillatory behavior near 7 = 1, 
as predicted in (81. 



k 


1 


4 


5 


10 


14 


21 


30 


55 


111 


140 


341 


7c,3,fc 


0.863 


0.852 


0.875 


0.851 


0.880 


0.857 


0.862 


0.860 


0.854 


0.891 


0.853 



Table 2. Values ^ c ,3,k of 7 at which some branches with k bound states cross 
the threshold R = 1. 



Contrary to the dimension 2, here some potentials with a higher number of bound states have 
a higher R than VCy d i- The corresponding curves in the 7 — R plane are flatter and flatter, and 
intersect at a sequence of increasing 7^. This sequence seems to accumulate at 1, in accordance 
to Helffer and Robert's result 18 , which predicts the existence of similar potentials with a R 
larger than 1 for every 7 < 1 in the nearly semiclassical regime. However, their study of an 
harmonic potential of varying width showed a highly oscillatory behavior of E(V) with respect 
to the width of the potential. Although our numerical methods do not permit us to investigate 
this regime (it would require a very large domain size, with a correspondingly large number of 
points, and a very poor conditioning of the eigenvalue problem), we expect the same behavior 
to occur. This means that the computation of the maximizing branches (here, those with 1, 5, 
14 and 140 bound states) becomes harder and harder. 

The relative energies of the branches display a greater variety than in the case d = 2, reflecting 
a more complicated energy landscape. Figure [8] shows the profiles of some critical potentials. 
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Note that the potentials in dashed lines can be regarded as "anomalous" versions of their 
similarly extended counterparts, and have a lower energy. As can be expected, the potential 
branches which are maximizers for some 7 have a profile that is decreasing in r. 






100 200 300 400 500 600 




1000 2000 3000 4000 5000 6000 



Figure 8. Shape of some critical points of E(V), at 7 = 1. The qualitative 
shape of the curves does not change much when varying 7. 



Table [3] displays the eigenvalue repartitions of a particular potential, the one with 140 bound 
states, at 7 = 0.88. Several patterns can be noted. First, for low values of k and I, the 
approximate relationship Xk+1,1 = A^+2 holds. This is because the associated eigenfunctions 
are localized close to 0. In this region, the potential can be approximated by a harmonic 
potential V(0) + ±V"(0)x 2 , which leads to the approximation X k>l = V(0) + y/W(0)(2k + l + §), 
explaining the relationship Xk+1,1 ~ Afc,z+2- This is only valid for small k and I, where the 
harmonic approximation is valid. When the eigenvalues are close to zero, another pattern 
emerges: the last negative eigenvalues have the same value of k + I, leading to a triangular 
pattern in Table [3| This seems to be true for the maximizing branches (1, 5, 14, 140), but not 
for the others (which tend to have a different ordering of the eigenvalues close to zero). We do 
not have an explanation for this. 
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A /- 7 


k = 1 


k = 2 


k = 3 


k = 4 


k = 5 


= 6 


k = 7 


k = 8 


I = 


-85 


-54 


-30 


-13 


-4 


-0.7 


-0.02 


+0.25 


I = 1 


-69 


-41 


-20 


-7.6 


-1.8 


-0.1 


+0.02 




I = 2 


-54 


-29 


-12 


-3.4 


-0.3 


+0.03 






I = 3 


-39 


-18 


-5.8 


-0.7 


+0.04 








I = 5 


-26 


-9 


-1.3 


+0.06 










I = 6 


-14 


-2.3 


+0.08 












I = 7 


-3.8 


+0.1 














jE 3. Eigenva. 


ue repartition : 


:or the 


branch with 140 bound states, 



7 



0.88, in units of 10 for readability. As in Section 2.2 Ay is the fc-th 



eigenvalue of equation ([5]). The positive eigenvalues are finite size effects and 
have no physical meaning, therefore we do not write them beyond the first one. 



Every potential we were able to compute was below the R = 1 threshold for 7 > 1. These 
results confirm the Lieb-Thirring conjecture that 7 Cj 3 = 1. However, the subtle behavior near 
7 = 1 means that there might exist maximizing potentials for 7 > 1 we were unable to find. 

4.4. The case d > 4. The results we obtained in the case d > 4 are similar to the case d = 3, 
with V-y^i as the maximizer for small 7, until it is outperformed by branches with a larger 
number of bound states, which all fall under the semiclassical limit before 7 = 1. However, 
the high cost of computation (the higher the dimension, the more ill-conditioned the eigenvalue 
problem is) prevents us from performing a systematic study. 

5. Conclusion 

In this paper, we used a maximization algorithm to investigate the best constants of the 
Lieb-Thirring inequalities, an important open problem of quantum mechanics. Discretizing 
the problem using finite elements, we were able to numerically compute critical points of the 
functional. Our main findings are the following: 

• In one dimension, the only critical point we found is the potential with exactly one 
bound-state V^y i 1. For all initial data, the algorithm seems to either converge to this 
potential, or to split in separating bumps. This supports the conjecture that these 
potentials are the only maximizers of the functional for ^ < 7 < |. 

• In two dimensions, the only critical points we found were radial. For all initial data, 
the algorithm seems to either converge to a radial potential, or to split in diverging 
bumps. Among the radial potentials we were able to compute, Ky,2,i was the maximizer 
for 7 < 1.16. After this point, the maximum corresponds to the semiclassical regime. 
The natural conjecture is that VC ^ i is the maximizer for 7 < 1.16, and therefore that 
7 Cj 2 « 1.16. 

• In three dimensions, branches of radial potentials with more than one bound state 
outperformed V^,d,i- This confirms the theoretical result from J8j, and provides new 
lower bounds (see Figure [7]). The natural conjecture is that there is a non-decreasing 
sequence of integers n 7 with ra 7 — > 00 as 7 — > 1 such that the maximizer of the functional 

) has n-y bound states. 

• In dimension d > 4, the results are similar to the three-dimensional case, and we con- 
jecture the same behavior. 

Beyond these results, the study hinted at a very rich and highly nonlinear behavior of the 
maximizers of Lieb-Thirring inequalities. This study is based on a simple numerical method (fi- 
nite element discretization). More involved computations could use a more appropriate Galerkin 
basis, such as the Gaussian basis functions used in quantum chemistry. This could allow for 
a more accurate computation of potentials with a larger number of bound states, and a more 
detailed exploration of the energy landscape. 
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On the theoretical side, the properties of the functional (P 7 dl remain unexplored. An open 
question is whether one can prove the existence of maximizers. The behavior of the maxi- 
mization algorithm is also an interesting question, with the separation in bumps particularly 
interesting. Finally, we believe that our method could be adapted to generalizations such as 
models with positive temperatures or positive density [7] , negative values of 7 [5] or polyhar- 
monic operators (—A)' + V (see [22]). 
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